{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests of variable cell minimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using JuLIP\n",
    "using JuLIP.Potentials: StillingerWeber\n",
    "using JuLIP.Solve: minimise!\n",
    "using JuLIP.Constraints: VariableCell\n",
    "using JuLIP.ASE: ASECalculator\n",
    "\n",
    "using PyCall\n",
    "using Plots\n",
    "using Polynomials\n",
    "\n",
    "@pyimport ase.units as units\n",
    "\n",
    "# @pyimport quippy.potential as quippy_potential\n",
    "# sw_pot = quippy_potential.Potential(\"IP SW\")\n",
    "# sw_calc_Q = ASECalculator(sw_pot)\n",
    "\n",
    "sw_calc_J = StillingerWeber()\n",
    "\n",
    "at = Atoms(\"Si\")\n",
    "#@assert energy(sw_calc_J, at) - energy(sw_calc_Q, at) < 1e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Speed test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Nsup = 10:20\n",
    "Natoms = 2*Nsup.^3\n",
    "t_J = []\n",
    "t_Q = []\n",
    "for N = Nsup\n",
    "    push!(t_J, @elapsed forces(sw_calc_J, at * (N, N, N)))\n",
    "    #push!(t_Q, @elapsed forces(sw_calc_Q, at * (N, N, N)))\n",
    "end\n",
    "\n",
    "plot(Natoms, [t_J], label=[:JuLIP :QUIP], marker=:o, \n",
    "    xlabel=\"N_atoms\", ylabel=\"Time / s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy volume plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E = Float64[]\n",
    "V = Float64[]\n",
    "P = linspace(-0.01, 0.01, 5)\n",
    "for p in P\n",
    "    at = Atoms(\"Si\")\n",
    "    set_calculator!(at, sw_calc_J)\n",
    "    set_constraint!(at, VariableCell(at, pressure=p))\n",
    "    minimise!(at, verbose=2)\n",
    "    V1 = det(defm(at))\n",
    "    @printf(\"p=%.3f, V1=%.3f\\n\", p, V1)\n",
    "    push!(E, energy(at))\n",
    "    push!(V, V1)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadratic fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p = polyfit(V, E, 2)\n",
    "V0, = roots(polyder(p))\n",
    "K = -p[1]\n",
    "@printf(\"quadratic fit: V0=%.3f, Bulk modulus K = %.3f GPa\", V0, K/units.GPa)\n",
    "\n",
    "scatter(V, E, marker=:o)\n",
    "plot!(v -> polyval(p, v), xlabel=:Volume, ylabel=:Energy, legend=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p = polyfit(V, collect(P), 1)\n",
    "V0, = roots(p)\n",
    "K = - V0 * polyval(polyder(p), V0)\n",
    "@printf(\"linear fit: V0=%.3f A^3, Bulk modulus K = %.3f GPa\", V0, K/units.GPa)\n",
    "\n",
    "scatter(V, P)\n",
    "plot!(x -> polyval(p, x), xlabel=:Volume, ylabel=:Pressure, legend=false)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
